home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / libry51.zip / LIBRY7A.DOC < prev    next >
Text File  |  1989-11-10  |  5KB  |  267 lines
  1. .de
  2. .pa
  3.          EXAMPLE PROGRAM ILLUSTRATING THE USE OF VECTOR INSTRUCTIONS
  4.  
  5.  
  6. $STORAGE:2
  7.       PROGRAM TEST
  8. C
  9. C  this program illustrates how to use the vector instructions
  10. C
  11.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  12.       PARAMETER (N=5)
  13.       DIMENSION AS(N,N),BS(N),CS(N),AV(N,N),BV(N),CV(N),JPIVOT(N)
  14. C
  15. C  create Vandermonde
  16. C
  17.       DO 100 I=1,N
  18.       BS(I)=FACT(I-1)
  19.       BV(I)=FACT(I-1)
  20.       DO 100 J=1,N
  21.       AS(I,J)=FLOAT(I)**(J-1)
  22.   100 AV(I,J)=FLOAT(I)**(J-1)
  23. C
  24. C  scalar solve
  25. C
  26.       CALL SCALAR(AS,BS,CS,JPIVOT,N,IER)
  27. C
  28. C  vector solve
  29. C
  30.       CALL VECTOR(AV,BV,CV,JPIVOT,N,IER)
  31. C
  32. C  list solution
  33. C
  34.       DO 200 I=1,N
  35.   200 WRITE(*,'(1X,I2,1P2E15.5)') I,CS(I),CV(I)
  36. C
  37.       STOP
  38.       END
  39.       SUBROUTINE SCALAR(A,B,X,JPIVOT,N,IER)
  40. C
  41. C  GAUSS ELIMINATION WITH FULL PIVOTING (SIMULTANEOUS EQUATIONS)
  42. C  SINGLE PRECISION VERSION WITHOUT VECTOR CALLS
  43. C
  44.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  45.       DIMENSION A(N,N),B(N),X(N),JPIVOT(N)
  46.       DATA SMALL/1.E-35/
  47. C
  48. C  CHECK FOR ERRORS
  49. C
  50.       IER=0
  51.       IF(N.LT.1) GO TO 900
  52. C
  53. C  N=1 SIMULTANEOUS EQUATIONS
  54. C
  55.       IF(N.GT.1) GO TO 100
  56.       A11=A(1,1)
  57.       IF(ABS(A11).LT.SMALL) GO TO 910
  58.       X(1)=B(1)/A11
  59.       GO TO 999
  60. C
  61. C  N>1 SIMULTANEOUS EQUATIONS
  62. C
  63.   100 N1=N-1
  64. C
  65. C  INITIALIZE THE PIVOT VECTOR
  66. C
  67.       DO 110 I=1,N
  68.   110 JPIVOT(I)=I
  69. C
  70. C  REDUCE THE MATRIX TO UPPER TRIANGULAR FORM
  71. C
  72.       DO 160 K=1,N1
  73.       K1=K+1
  74. C
  75. C  LOCATE THE LARGEST ELEMENT IN THE REDUCED MATRIX
  76. C
  77.       IP=K
  78.       JP=K
  79.       AMAX=ABS(A(IP,JP))
  80.       DO 120 I=K,N
  81.       DO 120 J=K,N
  82.       AIJ=ABS(A(I,J))
  83.       IF(AIJ.GT.AMAX) THEN
  84.         AMAX=AIJ
  85.         IP=I
  86.         JP=J
  87.       ENDIF
  88.   120 CONTINUE
  89.       IF(AMAX.LT.SMALL) GO TO 910
  90. C
  91. C  SWAP ROWS TO MOVE THE LARGEST ELEMENT INTO THE PIVOT POSITION
  92. C
  93.       IF(IP.EQ.K) GO TO 130
  94.       BK=B(K)
  95.       B(K)=B(IP)
  96.       B(IP)=BK
  97. C
  98.       DO 121 J=K,N
  99.       AKJ=A(K,J)
  100.       A(K,J)=A(IP,J)
  101.   121 A(IP,J)=AKJ
  102. C
  103. C  SWAP COLUMNS TO MOVE THE LARGEST ELEMENT INTO THE PIVOT POSITION
  104. C
  105.   130 IF(JP.EQ.K) GO TO 140
  106.       J=JPIVOT(JP)
  107.       JPIVOT(JP)=JPIVOT(K)
  108.       JPIVOT(K)=J
  109. C
  110.       DO 131 I=1,N
  111.       AIJ=A(I,K)
  112.       A(I,K)=A(I,JP)
  113.   131 A(I,JP)=AIJ
  114. C
  115. C  NORMALIZE THE ROW
  116. C
  117.   140 DO 150 I=K1,N
  118.       R=A(I,K)/A(K,K)
  119.       B(I)=B(I)-R*B(K)
  120.       DO 150 J=K1,N
  121.   150 A(I,J)=A(I,J)-R*A(K,J)
  122. C
  123.   160 CONTINUE
  124. C
  125. C  FINAL PIVOT ELEMENT
  126. C
  127.       ANN=A(N,N)
  128.       IF(ABS(ANN).LT.SMALL) GO TO 910
  129.       X(N)=B(N)/ANN
  130. C
  131. C  BACKSOLVE THE UPPER TRIANGULAR MATRIX
  132. C
  133.       DO 171 IN=2,N
  134.       I=N+1-IN
  135. C
  136.       BI=0.
  137.       DO 170 J=I+1,N
  138.   170 BI=BI+A(I,J)*X(J)
  139. C
  140.   171 X(I)=(B(I)-BI)/A(I,I)
  141. C
  142. C  PIVOT BACK TO THE ORIGINAL ORDER
  143. C
  144.       DO 180 I=1,N
  145.   180 B(I)=X(I)
  146. C
  147.       DO 181 I=1,N
  148.   181 X(JPIVOT(I))=B(I)
  149.       GO TO 999
  150. C
  151.   900 IER=1
  152.       GO TO 999
  153. C
  154.   910 IER=2
  155. C
  156.   999 RETURN
  157.       END
  158.       SUBROUTINE VECTOR(A,B,X,JPIVOT,N,IER)
  159. C
  160. C  GAUSS ELIMINATION WITH FULL PIVOTING (SIMULTANEOUS EQUATIONS)
  161. C  SINGLE PRECISION VERSION WITH VECTOR CALLS
  162. C
  163.       IMPLICIT INTEGER*2(I-N),REAL*4(A-H,O-Z)
  164.       DIMENSION A(N,N),B(N),X(N),JPIVOT(N)
  165.       DATA SMALL/1.E-35/
  166. C
  167. C  CHECK FOR ERRORS
  168. C
  169.       IER=0
  170.       IF(N.LT.1) GO TO 900
  171. C
  172. C  N=1 SIMULTANEOUS EQUATIONS
  173. C
  174.       IF(N.GT.1) GO TO 100
  175.       AA=A(1,1)
  176.       IF(ABS(AA).LT.SMALL) GO TO 910
  177.       X(1)=B(1)/AA
  178.       GO TO 999
  179. C
  180. C  N>1 SIMULTANEOUS EQUATIONS
  181. C
  182.   100 N1=N-1
  183. C
  184. C  INITIALIZE THE PIVOT VECTOR
  185. C
  186.       DO 110 I=1,N
  187.   110 JPIVOT(I)=I
  188. C
  189. C  TRANSPOSE MATRIX
  190. C
  191.       DO 120 I=1,N1
  192.       I1=I+1
  193.   120 CALL VSWP(A(I,I1),N,A(I1,I),1,N-I)
  194. C
  195. C  REDUCE THE MATRIX TO UPPER TRIANGULAR FORM
  196. C
  197.       DO 160 K=1,N1
  198.       K1=K+1
  199. C
  200. C  LOCATE THE LARGEST ELEMENT IN THE REDUCED MATRIX
  201. C
  202.       CALL VMAB(KK,A(1,K),1,(N-K+1)*N)
  203.       KK=KK+(K-1)*N
  204.       IP=(KK-1)/N+1
  205.       JP=KK-(IP-1)*N
  206. C
  207. C  SWAP ROWS TO MOVE THE LARGEST ELEMENT INTO THE PIVOT POSITION
  208. C
  209.       IF(IP.EQ.K) GO TO 130
  210.       BB=B(K)
  211.       B(K)=B(IP)
  212.       B(IP)=BB
  213. C
  214.       CALL VSWP(A(1,IP),1,A(1,K),1,N)
  215. C
  216. C  SWAP COLUMNS TO MOVE THE LARGEST ELEMENT INTO THE PIVOT POSITION
  217. C
  218.   130 IF(JP.EQ.K) GO TO 140
  219.       JJ=JPIVOT(JP)
  220.       JPIVOT(JP)=JPIVOT(K)
  221.       JPIVOT(K)=JJ
  222. C
  223.       CALL VSWP(A(JP,1),N,A(K,1),N,N)
  224. C
  225. C  NORMALIZE THE ROW
  226. C
  227.   140 DO 150 I=K1,N
  228.       R=A(K,I)/A(K,K)
  229.       IF(ABS(R).LT.SMALL) GO TO 150
  230.       B(I)=B(I)-R*B(K)
  231.       CALL VPIV(-R,A(K1,K),1,A(K1,I),1,A(K1,I),1,N-K)
  232.   150 A(K,I)=0.
  233. C
  234.   160 CONTINUE
  235. C
  236. C  FINAL PIVOT ELEMENT
  237. C
  238.       AA=A(N,N)
  239.       IF(ABS(AA).LT.SMALL) GO TO 910
  240. C
  241. C  BACKSOLVE THE UPPER TRIANGULAR MATRIX
  242. C
  243.       X(N)=B(N)/AA
  244. C
  245.       DO 170 IN=2,N
  246.       I=N+1-IN
  247.       I1=I+1
  248.       CALL VDOT(XI,A(I1,I),1,X(I1),1,N-I)
  249.   170 X(I)=(B(I)-XI)/A(I,I)
  250. C
  251. C  PIVOT BACK TO THE ORIGINAL ORDER
  252. C
  253.       CALL VMOV(X,1,B,1,N)
  254. C
  255.       DO 180 I=1,N
  256.   180 X(JPIVOT(I))=B(I)
  257.       GO TO 999
  258. C
  259.   900 IER=1
  260.       GO TO 999
  261. C
  262.   910 IER=2
  263. C
  264.   999 RETURN
  265.       END
  266. .ee
  267.